Application of Response Surface Methodology to Optimize Solid-Phase Extraction of Benzoic Acid and Sorbic Acid from Food Drinks

An experimental design was applied for the optimization of the extraction process of two preservatives, benzoic and sorbic acids (BA, SA), from food drinks. A simple, rapid, and reliable solid-phase extraction (SPE) method for the simultaneous extraction of these two preservatives and their determination by liquid chromatography with a diode array detector was considered. Box–Behnken design (BBD) was applied to both steps of the SPE process: (i) the sample percolation to ensure the retention of the totality of the acids by the silica-based C18 sorbent; (ii) the elution step to ensure desorption of the totality of the acids from the cartridge. Thus, the volume, pH, and flow rate of the sample, and the percentage of MeOH, volume, and flow rate of the elution solvent, were optimized. Sample volume and pH have a significant influence (p < 0.0001 and p = 0.0115) on the percolation yield. However, no effect was recorded for the flow rate (p > 0.05). Flow rate also has no significant effect on the elution efficiency. The proposed new solid-phase extraction method, which can be easily applied to routine monitoring of preservatives BA and SA in juice and soft drink samples, included 0.5 g of C18 sorbent, 1 mL of food drink adjusted to pH 1 and percolated at 4.5 mL min−1, and 1 mL of a solvent mixture composed of methanol/acidified water (pH = 2.6) (90:10, v/v) used in the elution step at a flow rate of 4.5 mL min−1. Validation of the SPE method and the technique of analysis were evaluated, namely, the accuracy, precision, detection, and quantification limits and linearity. Recovery percentages of benzoic and sorbic acids were above 95% with relative standard deviations lower than 1.78%. Detection and quantification limits were 0.177 and 0.592 µg mL−1, and 0.502 and 0.873 µg mL−1 for benzoic acid and sorbic acid respectively. Optimal conditions were applied to commercial fruit juices and soft drinks and a minimal matrix effect was observed. This method was compared with other SPE methods using oxidized activated carbon and multiwalled carbon nanotubes as adsorbents. The yields determined with these last two were low compared to those determined with our method.


Introduction
Food additives (colorants, preservatives, antioxidants, stabilizers, gelling agents, thickeners, flavors, flavor enhancers, sweeteners) are added to food to facilitate their manufacture, their conservation, and modify their nutritional qualities [1,2]. However, it becomes 1 mol L −1 ) were used. Certified standard solutions of BA and SA were obtained from Fluka (purity ≥ 99%). Standard stock solutions of 1 g L −1 of both acids were prepared in methanol/ultrapure water (40:60, v/v). Standard mixtures, ranging from 1 to 45 µg mL −1 , were prepared by appropriate dilution of the stock solution with acetate buffer/methanol (60: 40, v/v) [34]. The buffer solution (pH = 4.6) was prepared by mixing ammonium acetate and acetic acid in appropriate proportions. Silica-based sorbent with octadecyl functional group was acquired from Applied Separations Company, Allentown, PA, USA (3 mL; 500 mg, with a particle size of 40 µm and an average porosity of 60 A • ), and multiwalled carbon nanotubes (MWCNT) were purchased from Sigma-Aldrich (France). MWCNT were 5-9 µm in length and 110-170 nm external diameter. Activated carbon (AC) was purchased from Sigma (20-60 mesh, product number C3014).

Apparatus
The chromatographic analyses of BA and SA were carried out using a 1100 LC system (Agilent Technologies) equipped with a quaternary pump (model G1311A), degasser (model G7122A), automatic injector (model G1313A), diode array detector DAD (model G1315B), and Agilent ChemStation data processing software (Rev. B.02.01). Chromatographic separation was performed on an Agilent Zorbax SB-C8 analytical column (150 mm × 4.6 mm i.d., 5 µm particulate sizes), which were maintained at a temperature of 25 • C. The mobile phase consisted of 55.55% acetate buffer and 44.44% methanol adjusted to pH 4.5 with acetic acid [35]. The separation was performed in isocratic mode, the flow rate was 0.7 mL min −1 , and the injection volume was 10 µL. The identification was carried out using a DAD at a wavelength of 235 nm.

Solid-Phase Extraction Procedure
The solid-phase extraction of benzoic and sorbic acids was performed using a vacuum extraction module. Before loading samples, the cartridge was first activated with 10 mL of methanol and conditioned with 10 mL of UP water. After that, a volume of UP water spiked with BA and SA at 40 µg mL −1 was percolated through the cartridges. The volume, pH, and flow rate of the percolated sample were optimized to ensure the maximum retention of the two acids on the solid phase. The pH of samples was set between 1 and 3 with HCl or NaOH solutions. The solutes were then eluted with methanol/acidified water (1% acetic acid, pH = 2.6). The elution volume, percentage of MeOH, and the elution flow rate was also optimized to ensure desorption of all solutes from the SPE cartridge. These two SPE steps were optimized using the Box-Behnken design.

Experimental Design
Response surface methodology (RSM) was the technique applied for the optimization of the analytical procedure, and a Box-Behnken design was the model used [15,31,[34][35][36][37]. Three SPE factors were defined for the percolation and elution steps, and the number of experiments required was calculated by applying the following BBD equation: where N is the total number of experiments required; k is the number of factors (3 factors: X 1 , X 2 , X 3 ); and Cp is the center point. Accordingly, a design consisting of 15 experimental points including three center points was used to assess the effects of three variables and the interaction effects on responses by fitting the data to a polynomial model. The three most important parameters of each step (percolation and elution) were chosen as independent variables and named X 1 , X 2 , and X 3 , and set at three coded levels (−1, 0, +1). Design-Expert 13.0 software was used in the experimental design model building and data analysis.

Application of BBD for the Optimization of Percolation and Elution Steps
The two steps of the SPE method, sample percolation and acids elution, were optimized using the BBD. Specific volumes (factor X 1 ) of the standard solution (40 µg mL −1 ) were prepared and adjusted to the desired pH (factor X 2 ) and then percolated through the SPE cartridge at a flow rate (factor X 3 ) chosen by the percolation BBD. After passing the acids solution through the column, the recovered solution (reject) was analyzed by HPLC-DAD to determine the amount of acids adsorbed onto silica-based C18. Optimal conditions of this step were then applied when the next elution step was optimized. Indeed, the acids retained on the cartridge were eluted with methanol/acidified water at proportions (factor X1), volumes (factor X 2 ), and flow rate (factor X 3 ) experimentally chosen according to the BBD. The eluted solution (extract) was then analyzed by HPLC-DAD to determine concentrations of the eluted acids. These eluted acids give the extraction yield of the SPE method. Studied responses are the retention yields (Rr) and elution yields (Re) of benzoic and sorbic acids. These yields are calculated as follows: where Q 0 is the initial concentration of the acid in the percolated solution (40 µg mL −1 ); Q r is the concentration of the acid retained on the adsorbent during the percolation step; Q exp is the concentration of the acid determined in the recovered solution obtained after the percolation step; and Q e is the concentration of the acid determined in the extract after the elution step. These responses were predicted using the generalized second-order model given in the following equation: where R is the predicted response, β 0 is a constant, β i is the linear effect of variable X i , β ii is the quadratic effect of X i , β ij is the linear interaction effect of X ij , and ε is a statistical error term. Table 1 contains coded values and factor levels of the two BBDs applied to percolation and elution steps. The detailed matrix of BBD (for both percolation and elution steps) for the three-factor three-level design with the center is shown in Table 2.

Sample Preparation
The optimal conditions obtained with the optimized SPE method were then applied to commercial food drinks: fruit juices and soft drinks. Soft drinks were firstly degassed for 30 min and then filtered using 0.45 µm pore size membrane filters. While the samples of fruit juices are well homogenized to reduce to the maximum the pulp, if present. Then, 5 mL of each sample was diluted to 100 mL in ultrapure water and the pH was adjusted to 1 with hydrochloric acid. Finally, the SPE method was applied under its optimal conditions and the extract obtained was analyzed by HPLC-DAD. The same protocol described above for unspiked samples was applied to spiked samples with 1 mL of a standard solution containing BA and SA at concentrations of 10 µg mL −1 .

Comparative Study
The present SPE method was compared to other SPE protocols using oxidized activated carbon (AC) and multiwalled carbon nanotubes (MWCNT). The AC was oxidized with concentrated nitric acid HNO 3 at 25 • C by stirring for 24 h [36]. These two adsorbents were first oven-dried at 80 • C for 2 h and then packed in empty SPE cartridges (0.2 g, 3 mL, polypropylene) and to hold their packing sorbent in place; polypropylene upper and lower frits (20 µm porosity) were placed at each end of the cartridge. Carbon nanotubes were purified with an HCl solution (2 M) and washed with UP water. Finally, the cartridges are dried to remove any remaining water or solvents. Oxidized AC and MWCNT were used as SPE adsorbents using the following experimental conditions: conditioning with 10 mL MeOH and 10 mL ultrapure water; percolation of 1 mL of a standard solution of BA and SA (10 µg mL −1 ) adjusted to pH 1; and elution with MeOH/acidified water (90:10 v/v).

Effect of Experimental Parameters on Retention Yields
Parameters considered for the optimization of the percolation step are the sample volume (X 1 ), sample pH (X 2 ), and flow rate of the percolation (X 3 ). The experimental design and results are reported in Table 3. The experimental response studied is the retention efficiency, expressed in terms of percentage of retention (Rr%).

Model-Fitting and Statistical Analysis
The polynomial mathematical model developed for optimization is a second-degree model. Figures 1 and 2 showing the curves of predicted versus observed values confirm the goodness of fit (0.92 and 0.98). For each coefficient in the regression model, significance was assessed by the corresponding p values (p < 0.0026) [38]. The coefficients of determination (R 2 ) of the quadratic regression models vary between 0.98 ( Figure 1) and 0.99 (Figure 2), and values of predicted coefficients of determination vary from 0.92 to 0.96 for SA and BA, respectively, showing a reasonable agreement with the experimental results [39].

Model-Fitting and Statistical Analysis
The polynomial mathematical model developed for optimization is a second-degree model. Figures 1 and 2 showing the curves of predicted versus observed values confirm the goodness of fit (0.92 and 0.98). For each coefficient in the regression model, significance was assessed by the corresponding p values (p < 0.0026) [38]. The coefficients of determination (R 2 ) of the quadratic regression models vary between 0.98 ( Figure 1   Values of R 2 higher than 90% indicate a good fit between experimental values with those obtained by the models. The results obtained confirm that the actual values are very close to the predicted values (0.96 and 0.99). Mean values of 6.3 and 3.8 (p > 0.05) were obtained for the lack-of-fit test, suggesting that the model is reliable in predicting the re-  Values of R 2 higher than 90% indicate a good fit between experimental values with those obtained by the models. The results obtained confirm that the actual values are very close to the predicted values (0.96 and 0.99). Mean values of 6.3 and 3.8 (p > 0.05) were obtained for the lack-of-fit test, suggesting that the model is reliable in predicting the response and that it can predict variations within the response [37]. Figures 1 and 2 confirm that the curve of observed values as a function of predicted values perfectly fits the shape of a straight line; we note the close agreement that exists between the experimental results and the theoretical values predicted by the polynomial model [40][41][42][43][44].

Analysis of Significant Factors
The statistical significance of each term (linear, interaction, and quadratic) has been reported in Table 4 obtained from the analysis of variance. Each of the three terms can be considered for the second-order fit. The F-value of models of the first-order, two-way interactions, and pure quadratic, indicate that the models were significant at p < 0.05 for the percolation response of both acids. The lack-of-fit values, 3.8 and 6.3 associated with p-values of 0.21 and 0.14, were not significant due to the relative pure error. The results of the present study show that the model is significant with p < 0.05, and by this, the validity of the model is confirmed. Therefore, this model could work for percolation optimization. A summary of significant factors (p < 0.05) and their effect on the response variable are shown in Table 4. The results obtained show that the volume and pH have a significant influence (p < 0.0001 and p = 0.0115) on the retention yield. However, no effect was recorded for the flow rate (p > 0.05). The sample volume (X 1 ) played a primary role in improving percolation efficiency. Our results agree with other studies [38,42,45]. Overall, the statistical analysis suggests that the experimental values fit the models well, with good accuracy and reliability.

Mathematical Models
The coded values were analyzed using the Design-Expert software, with no transformation and quadratic model. The models showed adequate p-value and an insignificant lack of fit. The results obtained showed that the two models (R rBA and R rSA ) were statistically meaningful (their p-value < 0.000). Moreover, the coefficients R 2 for the responses R rBA and R rSA were, respectively, 0.99 and 0.98. These data reveal that the correlation was good, indicating a good fit of the models. The regression equations for the response variables in terms of coded factors are:

Effect of Interaction between Factors
The quadratic model obtained was used to calculate the response surface for each variable separately. Figures 3-5 show the response surface of Rr as a function of each pair of the independent variables. In Figure 3, the response model is mapped against sample volume (X 1 ) and pH (X 2 ), while the flow rate (X 3 ) is held constant at its central level. Examination of the results shows that retention yields increase when pH and sample volume decrease. As shown in Figure 3, the retention yields reached a maximum value at pH 1 and a sample volume of 1 mL.   Figure 4 shows the effect of sample volume (X 1 ) and flow rate (X 3 ) on the retention yields; pH is kept constant at 1. The diagnosis of response surfaces shows that the mutual interaction between the percolation volume and the flow rate is not significant (p > 0.001). The response surface plots show a significant effect of sample volume on percolation efficiency. However, no significant effect was observed for the sample flow rate (p > 0.001).
From response surface plots illustrated in Figure 5, we observe that the interaction between pH and the flow rate is not significant (p > 0.001).

Effect of Experimental Parameters on Elution Yield
After optimizing the percolation step, a second BBD with three factors and three center points was used to optimize the elution step. Independent variables were the percentage of MeOH % in the elution solvent (X 1 ), the elution volume (X 2 ), and the elution flow rate (X 3 ). Values of the elution yields of BA and SA did not follow the normal distribution. A logarithmic transformation of the results (log Re) was, therefore, necessary to apply the  Table 5.

Model Adjustment
The goodness of fit of the models was evaluated by the adjusted determination and the predicted determination coefficients. Polynomial mathematical models developed for optimization are second-degree models. The plot of predicted versus observed values (Figures 6 and 7) confirms the good fitting ability. Adjusted coefficients of determination of the models vary between 0.92 and 0.98, and values of predicted coefficients of determination vary from 0.82 to 0.92. Results obtained confirm that real values are very close to predicted values. Therefore, these models can be used to optimize responses due to the high correlation between observed and predicted values [37].   Table 6 reports the ANOVA results. The analysis of variance shows that the overall models are significant with p < 0.05, confirming the validity of these models for elution optimization. The results obtained show that the volume and % MeOH have significant effects (p < 0.05) on retention yield. On the other hand, no effect was recorded for the flow rate (p > 0.05). The p-values of linear coefficients (X1 and X2), interaction coefficients (X1X2), and quadratic coefficients (X2 2 ) were less than 0.05, which indicates significant effects on elution efficiency. Coefficients in the equations (R BA , R SA ) provide insight concerning the effects and the interaction between the factors that were studied (% MeOH, elution volume, and flow rate). The most important variables are the volume of the solvent, the percentage of methanol in the solvent, and their interactions. However, the elution flow rate has a negligible effect on the elution yield for both acids.   Table 6 reports the ANOVA results. The analysis of variance shows that the overall models are significant with p < 0.05, confirming the validity of these models for elution optimization. The results obtained show that the volume and % MeOH have significant effects (p < 0.05) on retention yield. On the other hand, no effect was recorded for the flow rate (p > 0.05). The p-values of linear coefficients (X1 and X2), interaction coefficients (X1X2), and quadratic coefficients (X2 2 ) were less than 0.05, which indicates significant effects on elution efficiency. Coefficients in the equations (R BA , R SA ) provide insight concerning the effects and the interaction between the factors that were studied (% MeOH, elution volume, and flow rate). The most important variables are the volume of the solvent, the percentage of methanol in the solvent, and their interactions. However, the elution flow rate has a negligible effect on the elution yield for both acids.  Table 6 reports the ANOVA results. The analysis of variance shows that the overall models are significant with p < 0.05, confirming the validity of these models for elution optimization. The results obtained show that the volume and % MeOH have significant effects (p < 0.05) on retention yield. On the other hand, no effect was recorded for the flow rate (p > 0.05). The p-values of linear coefficients (X 1 and X 2 ), interaction coefficients (X 1 X 2 ), and quadratic coefficients (X 2 2 ) were less than 0.05, which indicates significant effects on elution efficiency. Coefficients in the equations (R eBA , R eSA ) provide insight concerning the effects and the interaction between the factors that were studied (% MeOH, elution volume, and flow rate). The most important variables are the volume of the solvent, the percentage of methanol in the solvent, and their interactions. However, the elution flow rate has a negligible effect on the elution yield for both acids.

Mathematical Models
From the validation of the parameters studied (X 1 , X 2 , and X 3 ), the second-order polynomial was obtained, and it describes the response surface. The final equations, using the retention yield as a response for BA and SA, are, respectively: The validity of the model was determined by ANOVA. The coefficients of determination obtained were 0.92 and 0.98 for BA and SA, respectively, which indicates good agreement between the experimental and predicted values. The value of the F-test indicates that the second-order model was statistically significant (62.9 > 9.01).

Interaction Effects
The 3D response surface plots showing the elution yields against individual factors are shown in Figures 8-10, which illustrate the interaction of the volume of eluent with the percentage of MeOH, eluent flow rate with its volume, and the interaction of eluent flow rate with the percentage of MeOH, respectively. Results show an enhanced analytical response (Re %) when the percentage of MeOH is between 50 and 90 and the volume of eluent values is between 1 and 10 mL. Response surface plots reveal that the analytical response increased with decreasing volume of the eluent and increasing MeOH percentage.
Maximum elution yields were determined at eluent volumes and MeOH percentage close to 1 mL and 90%, respectively. Furthermore, the elution is enhanced when the elution flow rate is in the middle of its experimental values.

Determination of Optimal Conditions
Different factors can affect the SPE efficiency; therefore, their optimization through a multivariate approach is recommended, especially when these factors are correlated [46]. According to the results of the optimization of percolation and elution steps and desirability studies, experimental conditions given in Table 7 were chosen as optimal values for the SPE extraction of both studied acids. Similar results have been found in previous studies [41,47]. Table 7. Optimal conditions for the SPE of BA and SA.

Factor
Optimal Value R%

Method Validation
To consider all the concentrations of analyzed acids, during the optimization of the SPE method, several calibration ranges were considered (1-5, 1-15, and 25-45 µg mL −1 ). Three calibration curves were thus drawn for the two acids. Equations of the curves and the coefficients of determination R 2 , are given in Table 8. The R 2 values are greater than 0.9; we can therefore conclude that the regression models applied are linear. The LDD and LDQ were determined from regression lines determined for AB and AS (calibration range 1-5 µg mL −1 , repeated 3 times) and Y-intercepts were considered as blank responses. LDD and LDQ were calculated as shown below, and values are given in Table 8 [48].
where: k-factor of 3 and 10 for LDD and LDQ, respectively; s-standard deviation of the Y-intercept; and a-slope of the regression line. The precision of this method was studied [48]. For this purpose, five extractions of the ultrapure water spiked with the two acids at 5 µg mL −1 were performed using the optimum conditions described above. The repeatability of the SPE method was determined through the CV of the mean concentration obtained from the analysis of the five replicates achieved over a single day, while the reproducibility was calculated with results performed over different days and the CV of the mean concentration was calculated. The mean recovery percentages (~95%) showed a CV lower than 1.78% for the two acids, which highlights the high accuracy of the method. The data were also characterized by their relatively high precision. Indeed, CV values for intraday precision for AB and AS were lower than 1.48% and 2.34%, respectively. However, CV values for interday precisions for AB and AS were 1.24% and 1.85%, respectively.

Application to Real Samples
Robustness is the ability of an analytical method to provide small variations in results when subjected to controlled changes in application conditions (NF V 01-000). To verify the effect of a complex matrix and the robustness of our SPE method, all retained SPE conditions were assessed on real samples (three fruit juices and three soft drinks). Doped (with the two acids at 10 µg mL −1 ) and nondoped samples undergo the optimized SPE extraction method. The extracts obtained were analyzed by HPLC-DAD. The relative chromatograms are shown in Figures 11 and 12. As for extraction, yields are calculated as follows and shown in Table 9: where C 2 is the concentration of the acid in the extract obtained from a spiked sample, C 1 is the concentration of the acid in the extract obtained from a nonspiked sample, and C 0 is the concentration of the acid added to the sample (spiking level: 10 µg mL −1 ).

Application to Real Samples
Robustness is the ability of an analytical method to provide small variations in results when subjected to controlled changes in application conditions (NF V 01-000). To verify the effect of a complex matrix and the robustness of our SPE method, all retained SPE conditions were assessed on real samples (three fruit juices and three soft drinks). Doped (with the two acids at 10 µg mL −1 ) and nondoped samples undergo the optimized SPE extraction method. The extracts obtained were analyzed by HPLC-DAD. The relative chromatograms are shown in Figures 11 and 12. As for extraction, yields are calculated as follows and shown in Table 9: where C2 is the concentration of the acid in the extract obtained from a spiked sample, C1 is the concentration of the acid in the extract obtained from a nonspiked sample, and C0 is the concentration of the acid added to the sample (spiking level: 10 µ g mL −1 ).    The recoveries of benzoic and sorbic acids from real samples (Rs) are between 80 and 99.51% (Table 9). These yields correspond well to the model prediction and are very close to those determined during the SPE method optimization. Therefore, an insignificant effect of the matrix on the efficacy of our extraction method was observed. These results also show the importance of the application of an experimental design for the optimization of an experimental method or process. Figure 11. Chromatograms of the juice sample (a) and the juice sample spiked with 5 µ g mL −1 of BA and SA (b). Benzoic acid: 7.5 min; sorbic acid: 8.4 min.

Comparison Study
The present optimized SPE method was compared with an SPE method using an abundant and inexpensive adsorbent, oxidized activated carbon (AC), and with an SPE using multiwalled carbon nanotubes, a synthetic and expensive nanomaterial [47,48]. Table 10 shows the average recovery percentages with the relative standard deviations of the three SPE methods tested. Yields determined with both AC and carbon nanotubes adsorbents are low compared to those determined with silica-based C18. The oxidation of coal generally increases the carboxyl, lactone, phenolic, and carbonyl groups on the surface of this material. Carboxyl groups tend to form water clusters on the micropore openings in the AC surface, thus blocking benzoic and sorbic acids from entering the micropores and hindering the process of adsorption [39]. As for multiwalled carbon nanotubes, their morphology is nanotubular and multiwalled. This morphology is neither homogeneous nor regular, which does not ensure good reproducibility of extraction yields (RSD > 5%).

Conclusions
In the present study, a novel SPE method was developed for the enrichment of benzoic and sorbic acids from food drinks, to be analyzed by liquid chromatography. Therefore, this method could be easily adopted for routine monitoring of BA and SA preservatives in juice and soft drink samples. Compared to the usual method (International Standard ISO 22855-2008) [33], the proposed procedure is simple, fast, clean, and reliable. The extraction and preconcentration of these two acids before analysis increased the sensitivity of the detection method by reducing the matrix effect and concentrating these two acids in the final extract. On the other hand, the optimization of this SPE method by experimental design methodology (BBD) allowed for obtaining a robust, reliable, reproducible, and precise method with maximum extraction yield. The results showed good extraction yields (higher than 95%) under the following conditions: Conditioning of the C18 column was performed with 10 mL methanol followed by 10 mL of UP water. Then, 1 mL of sample, adjusted to pH 1, was percolated at a flow rate of 4.5 mL min −1 . Finally, the elution of acids was done with 1 mL of methanol/acidified water (pH = 2.6) (90:10, v/v) at a flow rate of 4.5 mL min −1 . Optimal conditions thus determined were successfully applied to commercial fruit juices and soft drinks. A slight matrix effect was observed since the calculated yields were close to those determined during the optimization. The present method was compared with SPE using an oxidized carbon AC and multiwalled carbon nanotubes MWCNT. Silica-based C18 presented better extraction yields. Thus, the proposed new solid-phase extraction method could be used for routine monitoring of preservatives BA and SA in juice and soft drink samples.